      program id5_beta
      implicit real*8 (a-h,o-z)
      parameter(ndec=12,npat=4096,ngrid=51)
c
      dimension pa(ndec,2),pe(ndec,2),tl(npat+1,3)
      dimension bgrid(ngrid),pz(ngrid),g(ngrid)
      dimension er1(ngrid,ngrid),er2(ngrid,ngrid)
      dimension in(npat,ndec)
      common /pat/in
c
      write(*,'('' ***** id5_beta *****'',20x,
     &/"    Beta=1 versus Beta<1"/)')
c
      open(unit=10,file='pat4096.data')
      do 10 n=1,npat
         read(10,*)idum,(in(n,k),k=1,ndec)
c         write(*,'(i5,x,12i2)')n,(in(n,k),k=1,ndec)
   10 continue
c
      do 15 j=1,ngrid
         bgrid(j) = 0.02*(j-1)
         if (j.eq.1) then                                
            pz(j) = 0.999d0
         else
            pz(j) = 0.5d0 + 0.01d0*(ngrid-j)
         endif
         g(j) = 2.d0*dlog(pz(j)/(1.d0-pz(j)))
   15 continue
c
      do 100 i=1,ngrid
       er2(i,ngrid) = 0.5d0              ! beta = 1
       gam = g(i)
       call fpx(gam,1.d0,pe)
       do 50 j=1,ngrid-1
         beta = bgrid(j)
         call fpx(gam,beta,pa)
         call fl(pa,pe,tl)
         call sort(npat+1,tl)
c
         fe = 1.d0
         fa = 0.d0
         do 30 n0 = 1,npat+1
            if (tl(n0,1) .ge. 0.d0) then
c                  write(*,'(2g12.4)')fe,fa
                  goto 31
            endif
            fe = fe - tl(n0,3)
            fa = fa + tl(n0,2)
   30    continue
   31    er1(i,j) = fa                    ! Prob(LL<x|Ha)
         er2(i,j) = fe                    ! Prob(LL>x|He)
   50  continue
       if (er1(i,j).gt.er2(i,j)) er2(i,j) = er1(i,j)
  100 continue
c
      write(*,'(//"Max of Type I and Type II Errors:"/)')
      write(*,'(7x,51f7.3)')(bgrid(j),j=1,ngrid)
      do 110 i=1,ngrid
         write(*,'(52f7.3)')pz(i),(er2(i,j),j=1,ngrid)
  110 continue
c
      end
c***********************************************************************
      subroutine fpx(gam,beta,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)

c      do 10 k=1,ndec
c         write(*,'(i2,2g12.4)')k,(px(k,j),j=1,2)
c   10 continue
c
      return
      end
c***********************************************************************
      subroutine fl(pa,pe,tl)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,npat=4096)
      dimension pa(ndec,2),pe(ndec,2),tl(npat+1,3)
      dimension in(npat,ndec)
      common /pat/in
c
      y1=0.d0
      y2=0.d0
      pt = 0.d0
      do 50 n=1,npat
         ta = 0.d0
         te = 0.d0
         do 10 k=1,ndec
            if (in(n,k) .eq. 1) then
               if (pa(k,1).gt.1.d-302) ta = ta + dlog(pa(k,1))
               if (pe(k,1).gt.1.d-302) te = te + dlog(pe(k,1))
            else
               if (pa(k,2).gt.1.d-302) ta = ta + dlog(pa(k,2))
               if (pe(k,2).gt.1.d-302) te = te + dlog(pe(k,2))
            endif
   10    continue
         tl(n,1) = ta-te
         tl(n,2) = dexp(ta)
         tl(n,3) = dexp(te)
         pt = pt + tl(n,3)
         y1 = y1 + ta*tl(n,2)
         y2 = y2 + te*tl(n,3)
c      write(*,'(i5,2x,5g12.4)')n,ta,te,tl(n,1),tl(n,2),tl(n,3)
   50 continue
      tl(npat+1,1) = 1.d2
      tl(npat+1,2) = 0.d0
      tl(npat+1,3) = 0.d0
c      write(*,'("pt = ",3g12.5)')pt,y1,y2
c      stop
c
      return
      end
c***********************************************************************
      subroutine sort(n,tl)
      implicit real*8 (a-h,o-z)
      dimension tl(n,3),a(n),b(n),c(n),indx(n)
c
      do 5 i=1,n
         a(i) = tl(i,1)
         b(i) = tl(i,2)
         c(i) = tl(i,3)
    5 continue
      call indexx(n,a,indx)
c
      do 10 i=1,n
         j = indx(i)
         tl(i,1) = a(j)
         tl(i,2) = b(j)
         tl(i,3) = c(j)
c      write(*,'(i5,x,3g12.4)')i,tl(i,1),tl(i,2),tl(i,3)
   10 continue
      return
      end

c***********************************************************************
      SUBROUTINE INDEXX(N,ARRIN,INDX)
      implicit real*8 (a-h,o-z)
      DIMENSION ARRIN(N),INDX(N)
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END


